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I,  GLOBAL  IONOSPHERIC  MODELING 


The  '/requirements  of  more  accurate  methods  for  assessing  HF  radio 
propagation  and  detection  conditions  has  intensified  th*'  efforts  to  produce 
better  models  of  the  ionosphere.  We  report  here  some  results  of  our 
continued  effort  in  this  area.  In  a previous  report  ^ we  described  some 
techniques  we  had  developed  and  some  computer  programs  written  toward 
the  construction  of  a three  dimensional  ionospheric  model.  This  model 
attempted  to  reduce  some  of  the  deficiencies  associated  with  models 
based  upon  our  current  knowledge  of  the  physical  and  chemical  processes 
occurring  in  the  ionosphere  by  the  inclusion  of  some  observations  of 
ionospheric  parameters. 

We  have  concentrated  our  effort  to  predicting  and  specifying  the  most 

(2) 

variable  ionospheric  region — the  F2  region.  Previous  work  has  indi- 
cated that  the  daily  variability  of  the  critical  frequency  of  the  F2  region, 
foF2,  is  two  to  three  times  that  of  the  critical  frequency  of  the  lower  iono- 
spheric regions.  The  techniques  we  employ  for  modeling  foF2  are  also  used 
in  modeling  hmF2,  the  maximum  height  of  the  F2  region  but  have  been  devel- 
oped only  from  studies  of  the  spatial  and  temporal  variability  of  foF2.  Studies 
using  ray  tracing  techniques  have  shown  the  inaccuracies  in  specifying  foF2 
along  a ray  path  are  responsible  in  most  part  for  Errors  in  predicting  HF 
propagation  conditions. 

It  is  shown  in  reference  1 that  our  approach  to  modeling  foF2  consists 
in  updating  or  modifying  a first  guess  ionosphere  with  observations  of 
foF2  both  at  irregular  locations  and  random  hours.  These  observations 
over  the  globe  are  extended  and  extrapolated  over  larger  regions  using  a ^ 

synoptic  mapping  technique  similar  to  that  used  by  meteorologists  for 
obtaining  weather  maps  when  only  sparse  observations  are  available.  The  * 

model  and  method  for  determining  the  vertical  electron  density  distribution 
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given  foF2  and  hn.F2  has  not  been  changed  and  will  not  be  described  here. 

We  have  used  several  sets  of  ionospheric  data  in  establishing  the  techniques 
and  in  verifying  the  accuracy  of  the  maps  that  were  generated. 

A.  The  1958  IGY  measurements  of  foF2  available  from  a network 

* of  157  vertical  incidence  ionosondcs  spread  over  the  entire 

globe. 

B.  Measurements  of  foF2  from  a set  of  36  vertical  incidence 
ionosrmdes  over  Europe  and  Asia  during  the  years  I960,  1962, 
1966  and  1968.  Data  was  also  available  from  the  same  network 
for  four  months  during  1964. 

These  data  span  a full  solar  cycle  and  allowed  one  to  test  the  modeling 
procedures  over  periods  of  both  high  and  low  solar  activity. 


V 


IONOSPHERIC  GRID  MODELING 


The  determination  of  foF2  on  a global  scale  consists  of  two  separate 
steps  which  can  be  modified  independent  of  each  other.  The  first  part  of 
the  modeling  procedure  consists  of  making  an  initial  estimate  of  foF2  at 
arbitrary  locations  and  times  over  a region  of  interest.  The  sparseness 
and  nonuniformity  of  the  available  data  requires  that  this  estimate  contain 
most  of  the  general  physical  features  of  the  ionosphere.  We  have  experi- 
mented with  several  first^guess  ionospheres  and  in  particular  problems 
have  used  more  them  one  in  producing  maps  for  series  of  hourly  analysis. 
The  most  general  accepted  predictions  of  foF2  are  the  monthly  iritf'  \an 

(3) 

maps  due  to  Jones  and  Gallet  . Their  model  describes  foF2  as  a 
spherical  harmonic  expansion  in  terms  of  0.  the  geographic  longitude,  X 
the  geographic  latitude  and  T.  the  universal  time.  The  diurnal  variation 
has  been  separated  from  the  spatial  variation  according  to 


foF2(0,  X,T) 


76 

= 1 


k = 0 


Dk(T)  Gk(0,X). 


(1) 


The  coefficients  of  the  expansion  have  been  related  to  the  12-month  running 
mean  Zurich  sunspot  number  so  that  the  model  predicts  median  ionospheric 
conditions.  The  monthly  median  model  is  used  whenever  no  other  analysis 
is  available.  When  ionospheric  maps  are  to  be  made  in  hourly  increments, 
we  have  used  as  the  first-guess  ionosphere  updated  maps  of  the  previous 
hour.  The  map  is  rotated  15°  of  geographic  longitude  to  the  west  as  a 
station  displaced  by  15°  to  the  west  would  be  expected  to  observe  a similar 
ionosphere  one  hour  later. 

Most  applications  of  our  predictions  require  that  the  ionosphere  be 
specified  on  a uniform  grid  coordinate  system.  The  output  of  our  analysis 
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is  always  given  in  geographic  coordinates  with  equal  spacings  in  latitude 
and  longitude.  We  have  found  that  a grid  spacing  of  15°  in  longitude  and 
10°  in  latitude  allows  us  to  preserve  most  of  the  important  ionospheric 
features  of  the  F2  region.  This  spacing  gives  us  a grid  which  is  25  X 19* 

In  some  areas,  most  notable  the  high  latitude  trough  and  the  region  of  the 
equitorial  anomaly  a smaller  spacing  in  latitude  may  be  necessary. 

In  our  original  work,  our  modeling  procedure- was  also  carried  out  in 
a uniform  grid  in  geographic  coordinates.  Since  many  features  in  the  F2 
region  tend  to  be  aligned  along  lines  of  constant  geomagnetic  latitude,  we 
have  chosen  values  of  the  grid  locations  (the  grid  is  still  25  X 19)  to  spaced 
15°  apart  in  geographic  longitude  and  10°  apart  in  geomagnetic  latitude. 
Although  this  is  not  an  orthogonal  coordinate  syp.tem,  it  has  no  effect  on  the 
modeling  procedure. 

The  first  guess  ionosphere  for  a given  universal  time  is  used  to  deter- 
mine the  value  of  foF2  at  each  of  the  grid  points  0 .,  X . and  at  the  location 

gi  mj 

of  each  observation  0^,  where  0^  is  geographic  longitude  and  is 
dipole  geomagnetic  latitude.  Using  this,  we  calculate  the  deviation  between 
first  guess  and  the  observation: 


= FG(0t  ^ T)  — OBS^ 


where  we  have  suppressed  the  indices  denoting  our  mixed  coordinate  system. 
We  also  calculate  for  the  N stations  reporting  the  average  residual 


Dk=*LD* 


and  the  average  square  residual  D_ 
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(4) 


N 

T i_  Z d) 

R = N 1=1 


The  influence  of  the  individual  deviations  D^is  extended  to  the  points  on  the 
grid  by  employing  a weighting  scheme.  The  corrections  to  the  first  guess 
field  at  a grid  point  for  an  individual  observation  can  be  written  as 

v 

V W<1U  w<2,ji  D*  <5> 

where  the  exact  form  of  the  W's  is  given  in  reference  (1).  The  weighting 
factors  depend  upon  north- south  and  east- west  geographic  separations  be- 
tween stations  and  grid  points.  The  total  correction  to  each  grid  point  is: 

(6) 


The  updated  grid  of  values  of  foF2,  Ujj  is  then  computed  using  these  correc- 
tions and  the  first  guess  field  using 


FG.. 


(7) 


The  updating  procedure  is  carried  out  a total  of  four  timus  with  the  updated 
field  U^j  becoming  a new  first  guess  field,  FG^  for  the  next  iteration. 

After  each  iteration,  a bilinear  interpolation  is  performed  between  the  four 
surrounding  g,  'd  points  of  each  observation  to  provide  a new  residual.  The 
magnitude  of  the  weighting  factors  is  reduced  by  a factor  of  two  after  each 
iteration,  but  is  not  allowed  to  decrease  below  a minimum  value,  If  we  define 
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* 


t 


t 


W 


(1) 


i-t 


j 


then  the  grid  correction  D-  can 


be  written  as 


L 

1=1 


-t=l  13 


l 


(9) 


In  order  to  examine  the  behavior  of  this  correction,  we  separate  out 
that  contribute  which  results  when  only  the  minimum  weight  Wo  is  used. 


M 

£ ’VD<-  + 

4=1  J 


N 

Wo  Z D, 

£=M+1 


L WiH  + (N — M)  W0 

i _ 1 ** 


(10) 


In  this  expression,  M observations  contribute  with  a weight  Wjj^>  WQ, 
and  N — M contribute  with  weight  Wo*  For  a given  local  time  and  season, 
W.j^is  dependent  upon  the  distance  between  grid  points  and  observations  and 
characteristically  goes  to  zero  between  2,500  and  7,500  kilometers  for  the 
first  iteration  and  between  300  and  900  kilometers  for  the  fourth  iteration. 
An  indication  of  the  size  of  the  correction  produced  by  Wo  not  being  zero  is 
given  in  regions  where  N=0,  that  is  no  nearby  observations.  Here  D.^  = 
the  average  residual  for  the  iteration.  In  almost  all  cases  | D_|  < 0. 1 Mhz . 
after  the  first  iteration.  In  addition  to  providing  an  average  correction  to 
the  entire  globe,  the  use  of  a non  zero  weight  insures  a smooth  behavior  of 
!>)..  at  the  boundary  of  regions  where  observations  do  not  exist.  If  N=1  and 
Wo=0  then  D..  = independent  of  separation  between  grid  point  and  obser- 

* J 
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vation.  This  condition  can  produce  a large  correction  for  distances  up  to 
7,500  kilometers. 

We  have  incorporated  into  the  analysis  program  the  capability  of  calcu- 
lating the  individual  deviations,  and  the  average,  D_u , and  rms  deviation 

R 

\j  for  groups  of  verifying  stations  not  used  to  update  the  first-guess  map. 

These  parameters  provided  an  indication  of  the  validity-  of  the  analysis  at 
points  away  from  the  observations  used  in  the  analysis.  We  present  here 
an  example  of  such  a study  in  which  the  157  stations  of  1958  were  randomly 
divided  into  two  sets  of  stations.  In  Table  X-l,  the  results  are  shown  for  four 
hours  during  January  24,  1958.  The  number  of  actual'  observations  available 
in  each  of  the  time  periods  is  comparable  for  both  station  groups.  Table  H-2a 
presents  a comparison  between  the  individual  observations  and  the  first  - 
guess  model  using  the  ITS  coefficients  and  the  final  updated  map  using  the 
modeling  stations  for  UT=1800  hours,  January  24,  1958.  Table  II-2B  presents 
the  same  information  for  the  verification  stations. 

We  can  see  from  Table  I that  when  observations  are  available  on  a 

global  basis  the  residuals  at  both  sets  of  stations  are  independent  of 

(4) 

time.  Studies  of  the  1958  data  have  shown  that  the  residuals,  and 


\j  D determined  from  the  modeling  stations  are  approximately  the  same  when 

the  number  of  stations  was  in  the  range  of  25-100.  It  can  be  stated  that  in 
regions  where  data  are  available,  the  modeling  procedure  yields  maps  that 
are  a good  representation  of  the  observations.  Figure  1-1  is  a contour  plot 
of  foF2  values  on  a cylindrical  equidistant  global  map  for  the  data  given  in 
Table  I-2&.  In  all  analysis,  the  general  shape  of  the  gradients  seen  on  the 
maps  that  used  daily  hourly  data  were  similar  to  the  monthly  median  for  the 
same  hour.  For  the  range  of  numbers  of  observations  used  in  the  analysis, 
the  general  shape  of  the  median  map  was  retained  and  the  values  of  map 
contours  were  raised  or  lowered  in  accordance  with  the  observations. 
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NUf3£*t  OF  STATiO-lS  IS  4 AVEF45E  RfcSIOUAL  F0F2  = _ -J,._5jL4  AVERAGE  RESIDUAL  F0F2  SQUARED,  = . 4.942 

(ITS  MODEL) 

NUM3ER  OF  3 T A T iO.N3  I.  4 AVERAGE  RE5ICJAL  F0F2  = .033  AVERAGE  RESIDUAL  F0F2  SQUARE*!  = .501 

(UPDATED  MAP  ANALYSIS) 
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Figure  1-1 


For  tiie  data  of  Table  1-2,  the  average  improvement  in  the  rms  residual 
for  the  sel  of  modeling  stations  is  1. 8Mhz.  and  for  the  set  of  verification 
stations  it  is  0. 7Mhz.  A more  detailed  study  of  the  improvement  in  the 
daily  specification  of  foF2  as  compared  with  the  median  ionospheric  structure 
ispresented  in  reference  4.  In  this  study,  the  improvement  resulting  from 
the  analysis  cf  daily  hourly  data  is  given  as  a function  of  various  latitudinal 
regions  and  seasons. 


IONOSPHERIC  PREDICTIONS 


tx * 


The  first  guess  foF2  field  each  hour  can  be  considered  the  prediction 
for  that  hour.  The  advantage  of  the  ITS  monthly  median  modei  is  that  it 
can  be  used  as  a prediction  for  arbitrary  time  periods.  The  chief  dis- 
advantage in  using  the  median  prediction  is  that  it  does  not  reflect  the 
day-to-day  variability  of  the  F2  region.  The  updated  maps  based  upon 
our  modeling  procedure  represent  a more  accurate  description  of  the  F2 
layer  for  an  hour  in  which  observations  of  foF2  are  available.  We  have 
studied  two  methods  by  which  these  analysis  of  past  data  may  be  used  to 
make  a prediction  for  future  ionospheric  behavior. 

The  first  method,  which  we  mentioned  when  describing  the  first  guess 
ionosphere  is  to  use  the  updated  analysis  made  one  hour  in  advance  of  the 
desired  prediction.  The  prediction  for  foF2  is  obtained  from  the  grid 
values  by  using. 

f°F2(0gi.  *mj.  T+I)  = foF2(0g.+ j,  Xm.,T)  (11) 

where  the  variables  have  already  been  described.  Since  the  modeling 
procedure  is  based  upon  the  calculation  of  angular  separations  between 
grid  points  and  observation  locations  in  an  orthogonal  coordinate  system, 
we  must  perform  transformations  between  this  mixed  coordinate  system 
and  geographic  coordinates.  Both  coordinate  systems  are  shown  in  the 
following  figure. 
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IForld  map  of  magnetic  dipolo  lalUudt. 


Figure  1-2 


The  grid  used  in  expression  (11)  is  specified  by  equal  spacing  in  magnetic 

dipole  latitude,  To  carry  out  all  transformations,  we  store  the  grid 

of  values  of  geographic  latitude  corresponding  to  the  regularly  spaced 

coordinates,  0 .,  X ..  That  is 
gv  mj 


i 


X ..=  H(0  X •) 
gij  vvgi'  my 


(12) 


new 


This  grid  is  determined  using  the  equations  which  specify  the  transformation 
of  one  spherical  coordinate  system  to  another  by  a rotation  matrix.  In  our 
case,  we  use  the  equation 


sin  = (R**  cos  0 . + R,o  sin0  .)  cos 

m]  gl  gl  g 


(13) 


where  ^13  = COfl  2°  cos  291.0° 

R23  = C0S  78,  3°  Sin  291,0° 
R33  = sin  78.  3° 


I 


For  each  grid  value  0g.,  we  choose  one  of  the  two  solutions  of  equation 
(13)  for  X ...  As  can  be  seen  from  the  world  map  of  figure  2,  not  all  values  of 
0 ^ . lead  to  a corresponding  Agij.  The  results  for  the  25  X 19  grid  of 

gi  mj 

equal  spacing  of  0 , is  given  in  Table  1-3,  where  a value  of  999.  0 indi- 
cates that  a coordinate  does  not  yield  a transformation.  The  modeling  pro- 
cedure is  carried  out  in  the  regular  mixed  coordinate  system  0gi,  Xmj  with 

the  weights  of  observations  at  the  grid  points  given  in  terms  of  the  angular 

obs  obs 

separations  (0g^  — 0g£  ) and  (^gij  — )•  In  order  to  determine  foF2  at 

arbitrary  locations,  a linear  interpolation  in  magnetic  dipole  latitude  is 
performed  on  the  great  circle  of  constant  geographic  latitude.  These  inter- 
polations are  carried  out  after  each  iteration  step  in  order  to  calculate  the 
new  residuals  at  station  locations  and  at  the  end  of  the  analysis  to  determine 
a final  grid  map  in  geographic  coordinates. 


Several  tests  were  performed  to  determine  the  differences  between  this 
modeling  procedure  and  one  in  which  only  geographic  coordinates  are  used. 
If  no  observations  are  used,  the  only  possible  change  in  the  final  map  must 
be  due  to  the  successive  interpolations  performed  on  the  median  ionosphere 
between  the  two  coordinate  systems.  This  was  observed  to  be  of  the  order 
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of  a few  tenths  of  megahertz  at  the  equator  where  the  north-  south  gradients 
of  foF2  are  the  largest.  When  a large  sample  of  observations  was  used  in 
the  analysis,  the  final  maps  produced  using  the  two  different  coordiante 
systems  were  in  agreement  with  each  other -to  better  than  5%  except  at  the 
equator  were  the  differences  between  grid  points  were  less  than  10%. 

The  second  method  for  making  a predictioivfor  the  behavior  of  F2  layer 
.makes  use  of  both  the  immediate  past  observations  of  foF2  and  the  modeling 
procedure  we  use  for  updating  current  maps;  We  assume  that  data  were 
available  at  particular  locations  at  the  same  hour  for  which  we  wish  to  make 
a prediction  for  the  preceding  five  days.  We  use  the  short  term  prediction 
scheme  for  foF2  described  by  Rush  and  Gibbs  The  five  day  weighted 
mean  predictica,  Do,  is  given  by 

D0  = (5D_1  + 4 D_2  + 3 D-3  + 2 D_4  + D_5)_l  (14) 

15 

where  Dj  is  the  value  of  foF2  for  the  i^  day  preceding  the  prediction  day. 

All  the  Dj's  are  for  the  same  local  time  at  a station  location.  The  set  of  E^'s 
represent  a prediction  at  particular  locations  on  the  globe.  We  use  our 
modeling  procedure  to  extend  the  influence  of  this  set  of  isolated  predictions 
to  other  points  on  the  globe.  The  ITS  median:  ionosphere  is  updated  by  the 
individual  predictions  to  produce  a global  prediction  grid. 

The  two  prediction  techniques  have  been  tested  using  foF2  data  available 
for  the  years  I960,  1962,  1 966  and  1968.  The  stations  used  for  I960  are 
shown  on  the  global  map  of  figure  1-3  and  listed  in  table  1-4.  The 
results  we  present  here  are  for  an  entirely  different  method  of  selecting  the 
sets  of  stations  to  be  used  for  modeling  and  verification.  With  the  1958  data, 
we  choose  these  sets  in  a random  manner  so  that  each  set  of  stations  would 
give  approximate  uniform  coverage  of  the  globe.  For  thio  selection  pro- 
cedure, we  observed  that  the  updating  analysis  produced  an  foF2  which 
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Figure  1-3 


’Verification  Stations 


Modeling  Stations 


Station 

Code 

Longitude 

Latitude 

041 

12.500 

41.800 

048 

7.600 

48.100 

050 

10.100 

51.600 

051 

-.600 

51.500 

052 

14.600 

50.000 

05 5 

13.400 

54.600 

138 

23.600 

38.000 

142 

44.800 

41.700 

145 

20.500 

44.800 

149 

39.700 

47< 200 

152 

21.200 

52.200 

155 

37.300 

55.500 

156 

44.300 

56.100 

158 

17.60G 

59.800 

159 

24.60 J 

60.500 

160 

30.700 

60.000 

164 

18.800 

64.700 

166 

26.600 

67.400 

167 

20.400 

67.800 

168 

33.000 

69.000 

236 

51.400 

35.700 

237 

58.300 

37.900 

256 

61.000 

56.700 

26t> 

66.500 

66.500 

343 

76.900 

43.200 

356 

84.900 

56.500 

Station 


Code 

Longitude 

Latitude 

352 

104.000 

52.500 

424 

121.200 

25.000 

426 

127.800 

26.3Q0 

431 

130.600 

31.200 

462 

129.700 

62.000 

535 

139.500 

35.700 

539 

140.100 

39.700 

545 

141.700 

45.400 

547 

143.000 

47.000 

664 

-173.400 

64.400 

TABLE  1-4 
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significantly  reduced  the  residuals  at  each  of  the  two  sets  of  stations.  For 
the  I960  data,  we  chose  as  our  modeling  stations  all  those  that  were  west 
of  90°E  longitude  and  as  our  verification  stations  all  those  that  were  east 
of  90°E  longitude.  This  selection  should  give  some  indication  of  how  good 
the  modeling  procedure  is  in  regions  far  removed  from  observations. 

The  first  step  when  using  the  five  day  weighted  mean  (FDWM)  scheme 
is  to  produce  a map  using  the  median  model  as  a first  guess.  We  have  used 
the  data  of  the  24  hours  of  January  14,  I960  as  a test  case.  For  this  test 
period,  approximately  19  stations  reported  FDWM  data  for  modeling  (M)  and 
10  stations  reported  FDWM  data  for  verification(V).  The  results  of  this 
analysis  are  presented  in  Table  1-5.  We  first  note  that  the  values  of  the  RMS 
residuals  for  the  first  guess  for  all  of  the  hours  for  both  sets  of  stations  is 
less  than  that  for  daily  data.  Here  the  average  RMS  is  approximately  1 Mhz 
whereas  for  the  daily  data  of  1958  it  was  greater  than  2 Mhz.  The  FDWM 
data  is  in  much  better  agreement  with  the  median  model  than  the  daily  data. 
After  the  updating  procedure  was  applied,  the  RMS  residual  for  the  modeling 
stations  was  reduced  to  less  than  0.5  Mhz  while  the  RMS  residual  for  the 
verification  stations  shows  a smaller  reduction  for  approximately  70%  of  the 
hours  studied.  This  improvement  is  realized  inspite  of  the  large  north- south 
gradients  (a  10  Mhz  change  in  40°  latitude)  observed  in  the  verification  data 
including  stations  reporting  from  within  the  equitorial  anomaly  region. 

Table  1-6  shows  the  results  for  one  of  the  hours  for  each  station.  This  indi- 
cates that  it  is  possible  to  obtain  an  improvement  in  our  knowledge  of  foF2 
at  a particular  time  in  areas  where  observations  are  not  readily  available. 
The  final  grid  maps  for  each  of  the  hours  can  now  be  used  as  a prediction 
field  for  the  following  day. 

We  use  the  hourly  measurements  of  foF2  for  January  15,  i960  for  both 
modeling  and  verification  stations  to  study  the  three  methods  of  making  a 
prediction  for  foF2.  The  two  grid  field  methods  we  have  just  described 
can  each  be  used  in  several  ways.  The  updated  FDWM  field  could  be  used 
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.RMS  Residuals 

Five  Day  Weighted  Mean  Data 
January  14,  i960 


Universal  Time 

First  Guess  (Median  Model) 

Final  Map 

M 

V 

M 

V 

000 

1.0 

1.4 

0.3 

1.7 

1.1 

1.2 

0.5 

1.1 

200 

1.1 

0.7 

0.5 

0.8 

1.0 

0.9 

0.5 

1. 1 

400 

1.0 

0.8 

0.6 

0.9 

0.9 

0.7 

0.5 

0.7 

600 

0.9 

0.8 

0.4 

0.7 

0.8 

1.0 

0.5 

0.9 

800 

0.8 

1.0 

0.4 

0.8 

0.7 

0.9 

0.4 

0,9 

1000 

0.9 

0.7 

0.3 

CM 

• 

1.0 

0.8 

0.5 

0.9 

1200 

1.0 

0.9 

0.5 

0.8 

1.0 

0.9 

0.4 

0.7 

1400 

0.9 

0.  8 

0.5 

0.6 

0.9 

0.9 

0,3 

0.8 

1600 

0.9 

0.  8 

0.5 

0.7 

0.8 

0.8 

0.3 

0.6 

1800 

0.9 

0.9 

0,4 

0.7 

0.9 

0.  8 

0.  3 

0.6 

2000 

0.9 

0.9 

0.6 

0.  0 

1.0 

1.0 

0.2 

1.2 

2200 

0.9 

1.0 

0.5 

1.0 

1.1 

1.2 

0.4 

1.3 

Table  1-5 
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as  a prediction  for  more  than  one  day  in  the  fir.  \re  and  the  maps  obtained 
by  modeling  in  the  mixed  coordinate  system  could  be  rotated  nl5°  for  a n**1 
hour  in  the  future  prediction.  We  present  here  only  the  closest  time  internal 
application  of  both  schemes.  In  addition,  we  present  the  prediction  obtained 
using  the  median  model  for  comparison.  This  model  requires  the  smoothed 
sunspot  number  to  calculate  a value  of  foF2.  The  two  grid  field  maps  are 
used  without  any  modification  based  upon  any  additional  knowledge  of  iono- 
spheric activity. 

The  RMS  residuals  for  the  three  methods  for  the  24-hour  period  are 
given  in  Table  1-7.  As  a prediction,  there  should  be  no  difference  between  the 
modeling  and  verification  stations  using  the  median  model.  The  consistently 
large  residuals  from  the  verification  stations  must  be  due  to  local  dynamic 
effects  in  the  ionosphere.  For  this  reason*  we  only  compare  the  predictions 
for  each  set  of  stations  separately.  For  almost  all  of  the  hours,  there  is  no 
significant  difference  between  the  three  predictions  for  the  verification  stations,. 
For  the  modeling  stations,  there  is  some  slight  improvement  over  the  median 
model  by  the  two  grid  predictors  for  a few  of  the  hours.  But  even  here,  no 
systematic  behavior  is  observed.  Other  modifications  to  the  modeling  pro- 
cedure might  improve  the  prediction  capability  of  the  five-day  weighted 
mean  data.  A more  detailed  study  of  the  prediction  capability  of  our  modeling 
procedure  is  given  in  reference  4.  Here,  it  is  shown  that  one  can  improve 
upon  the  median  specification  of  the  F2  region  in  the  midlatitude  ionosphere. 

As  we  stated  earlier,  the  analysis  of  the  F2  region  has  been  coupled  with  a 
vertical  electron  density  model  to  produce  a complete  three  dimensional 
ionospheric  model  which  is  currently  being  used  in  a variety  of  communi- 
cation applications. 


RMS  Residuals  for  Various  Predictions  of  foF2 


January  15  January  15  January  15 

First  Guess  First  Guess  First  Guess 


Universal  Time 

FDWM 

ITS 

Rotated  1 Hour 

M 

V 

M 

V 

M. 

V 

000 

1.1 

2.3 

1.5 

2.4 

1.5 

2.3 

1.4 

1.8 

1.6 

2.0 

0.8 

1.7 

200 

1.0 

1.6 

1.0 

1.4 

1.1 

1.3 

1. 1 

1.7 

1.2 

1.7 

1.1 

1.7 

400 

1.0 

1.9 

1.2 

1.7 

1.0 

2.1 

0.9 

1.6 

1.4 

1.6 

2.0 

2.5 

600 

1.4 

2.4 

1.7 

2.3 

2.6 

2.2 

1.5 

1.9 

1.6 

1.7 

1.3 

1.8 

800 

1.4 

1.8 

1.4 

1.0 

0.8 

2.2 

l.l 

1.7 

1.1 

1.1 

0.7 

2.0 

1000 

1.0 

1.3 

1.0 

1.2 

1.1 

1.5 

0.9 

1.3 

1.3 

1.6 

1.2 

1.4 

1200 

0.7 

2.6 

1.1 

2.6 

0.8 

2.4’ 

0;6 

2.2 

1.1 

2.3 

0.9 

2.1 

1400 

1.5 

2.5 

1.8 

2.6 

1.5 

2.1 

0.9 

2.0 

1.2 

2.1 

0.8 

1.7 

1600 

0.7 

2.7 

0.8 

2.8 

0.9 

2.5 

0.7 

1.5 

0.8 

1.6 

0.6 

1.6 

1800 

1. 1 

1. 1 

0.7 

1.2 

0.5 

1.1 

0.7 

1.1 

0.6 

1.1 

0.7 

1.1 

2000 

0.5 

1.1 

0.7 

1.2 

0.6 

1.1 

0.8 

1.2 

0.9 

1.1 

0.9 

1.1 

2200 

1.0 

1. 1 

0.9 

1.1 

0.9 

1.4 

1.3 

1.6 

0.9 

1.5 

Table  1-7 
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H.  STUDIES  IN  POLAR  IONOSPHERIC  MODELING 


In  this  section  we  will  describe  some  attempts  at  improving  our 
knowledge  of  the  F2  region  of  the  polar  ionosphere  through  the  use  of 
empirical  models.  Reference  (1)  describes  some  recent  work  per- 
formed toward  this  end  and  our  approach  reflects  the  same  requirements 
that  are  described  in  that  report.  Our  prime  requirement  is  the  accuracy 
of  specifying  the  behavior  of  foF2  in  the  polar  region  as  a function  of  time 
and  other  physical  parameters.  This  model  can  then  be  used  with  other 
ionospheric  parameters  to  predict  the  three  dimensional  distribution  of 
electron  density  and  aid  in  the  evaluation  of  systems  employing  radio 
propagation  through  the  ionosphere.  The  work  described  here  has  concen- 
trated on  improving  our  knowledge  of  foF2  and  has  made  use  of  severa. 
of  the  presently  used  empirical  models  of  the  ionosphere^^. 

The  most  widely  used  of  these  models  is  the  ITS  monthly  median 
(2) 

model'  which  represents  the  median  behavior  of  foF2  (and  other  ionospheric 
parameters)  in  terms  of  a spherical  harmonic  expansion  in  terms  of  uni- 
versal time  and  spatial  coordinates.  The  representation  is  of  the  form 

N 

f°F2  * ^ Dfc  (T)  Gfc  (0,  X)  where  0 and 
k=0 

\ are  the  geographic  longitude  and  latitude  and  T is  the  universal  time. 

Here  the  diurnal  variation  has  been  separated  from  the  spatial  variation 
and  the  cutoff  parameter  of  the  series,  N is  equal  to  75.  The  time  dependent 
functions  are  expandedin  a Fourier  series  containing  six  harmonics: 

G 

D^T)  = + I (Aj^  cos  jT+  sin  j'J). 

j=l 
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This  results  in  a function 


N 

foF2  = J Ao(k)Gk(0,X)  + 

k=0 

N 6 

I I (A,(k)cos(jT)Gk(0,X)+  B,(k)sin(jT)Gk(#,X)] 
k=0  j=l  J J 

where  the  G(0,  X)  are  linear  combinations  of  surface  spherical  harmonics. 
The  explicit  arguments  of  the  G functions  were  chosen  to  be  the  modified 
magnetic  dip  angle  which  is  a function,  of  latitude  and  the  magnetic  dip 
angle  I,  given  by 


I = tan"1 


where 


Bx,  By,  Bz  are  the  north,  east  and  vertical  components  of  the  earth 
magnetic  field  vector.  The  modified  magnetic  dip  angle  is  then 


The  988  coefficients  were  determined  using  the  hourly  observations  of 
foF2  that  were  available  from  a network  of  vertical  incidence  ionosondes 
for  tlie  IGY  1958.  The  use  of  additional  data  from  1964  allov/ed  for  the 
determination  of  the  mean  sunspot  number  dependence  of  the  coefficients. 

The  methods  employed  in  the  development  of  the  monthly  median 
model  and  the  distribution  of  data  tended  to  give  more  weight  to  the  low 
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latitude  representation  of  the  ionosphere  that  than  to  the  high  latitude.  It 

(3)  < 

has  been  shown  that  the  true  physical  variations  about  the  median  values 

of  the  F2  data  in  the  polar  region  is  much  greater  than  in  the  midlatitude 

(3) 

regions;  This  has  resulted'  in  an  improved  model  for  representation  of 
foF2  in  the  polar  region.  This  model  describes  the  behavior  cf  foF2  in 
terms  of  the  same  coordinates  and  parameters  and  does  not  contain  a pro- 
vision for  obtaining  a piecewise  continuous  value  of  foF2  oyer  the  entire 
globe. 

A more  recent  empirical  model  of  the  polar  ionosphere  is  the  AFCRL 
model  of  Elkins  and  jRush^.  This  predictive  model  of  t'ue  ionosphere  is 
based  upon  various  kinds  of  data  including  1958  and  1964  vertical  incidence 
ionosonde  measurements,  and  optical  and  satellite  observations.  The  repre- 
sentation of  the  different  regions  of  the  ionosphere  range  from  predominantly 
statistical  models  to  purely  deterministic  models.  Testing  of  the  entire  model 
has  revealed  that  the  greatest  need  for  improvement  is  in  the  specification  of 
foF2.  We  summarize  the  model  for  foF2  in  the  following  pages. 
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The  variables  and  physical  parameters  used  to  define  foF2  are: 

t = local  time  of  day  in  hours 
tjH  s corrected  magnetic  time  in  hours 
D = day  of  the  year  (i.  e. , 1-366) 

X = geographic  latitude  (north) 

<P  = geographic  longitude  (east) 

Xm  - corrected  geomagnetic  latitude  (north) 
s = sunspot  number 
kp  = planetary  magnetic  activity  index 
0)  a 2TT/24 
u = 2tt/365 

The  value  of  foF2  is  given  by 


f0F2  = B(D,  s)  + I(s) 

C1  + c2 

r "~ 

3 

A0  + 2 ) Ancoslnuitt+^J) 

n=l 

where 


B(D,  s)  = C3N(s)cos(u(D-J  p)) 

N(s)  = +1  si  100 
N(s)  = -1  s < 100. 

. 2 
I(s)  = c4+  c5s+  c^s 
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agH(spsD) 


m(kp)  = 1 

m(kp)  = 2 

m(kp)  = 3 
m(kp)  = 4 
m(kp)  = 5 
m(kp)  = 6 
m(kp)  = 7 


kp<  0.3 

0.  3 s kp  <1.3 

1.  3 <kp  < 2.3 

2.  3 < kp  < 3.3 

3.  3 <kp  <4.3 

4.  3 s:  kp  < 6.  3 
6.  3 < kp 


G(\n,  m(kp),  D)  = Pl  + p2Xfn(3X  + X^fX^  l.  5S) ) 


X = Xm°  (min)Xm°  (max) 

^ A 

S = Xm  (min)  + Xm°  (max) 

\n° (ltlax>  c % * <*5  c°*  (2uD) 

Aq 

\n  (min)  = ql  + q2  m(V  + q3  cos  (^(D+  e,)) 
P1  “ ri  + r2  co8(2u{D'+-02)) 
p2  a,r3  + r4  cos(2u(D+02)) 


H(CPf  D)SQ 


H(CP,D)  = 


l 


y_  d 

n 


n 


nsO 


165*  <cp  195*  (Pacific  Zone) 
195  «PsS345?  (American  Zone) 
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-15°  < «p  165°  (European  Zone) 


4 

H(CP,D)  = Y 0JnDn 
n=0 

Ax  = A24  for  X * 58° 

A1  = A24(90‘  " X)/32‘  for  X > 58° 

3. 

A24  = bo  + 2 I-  ncos(nu(D  + a) ) 
n=l 

3 

A2  = do  + 2 Z cos(nu(D+  6n) ) 
n=l 

3 

A3  = eo  + 2 Z «n  cos(nu(D  + en)) 
n=l 

3 

*1  = fo  + 2 7 fn  cos(nu(D  + pn)) 
n=l 

3 

*2  = go  + 2 Z gn  cos(nu(D+  \) ) 
n=l 

3 

i3  = hQ  + 2 £ hn  cos(nu(D+  Yn)) 
n=l 

Auroral  Oval  Correction  to  foF2 

for  6 * t < 18  (magnetic  day) 

A)2 

f0F2«  = f0F2(l-C)e^~ ~ 

where  C = 1 - .35  cos(u)(t  - $j))  (l  + cos(u(D  + 8)))/2 

A = 80-1.  2k^  (magnetic  latitude  of  auroral  oval 
p center  line  for  day,  Degs.) 


for  -6  ^ tm  < 6 

f0F2'  = f0F2+[kp/4]cos(UJtm)e 


-(Wei2 

6 


where  {3  = 72  - 1,  8kp  + 5.1  cos  (U)^-  1) ) 


(magnetic  latitude  of 
auroral  oval  center 
line  for  night,  Degs. ) 
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Constants  for  AFCRL  foF2  Polar  Model 

Cj  = 4.8 

3Lq  = 0#  0 

qj  = 88. 5 

vQ  = 0.521 

= 0.42 

aj  = -0.439xl0_1 

q2  = ~2,  5 

Vj  = -0.  347x10* 

C3  = 0.  6 

a2  = 0.386xl0*2 

q3  = 5.0 

v2  = 0.  316x10* 3 

c„  = 1.0 
4 

a3  = -0.424 

q4  = 55.0 

v3  = -0. 133x10* 

v.  = 0.142xl0'8 
4 , 

C5  a 0.12 

a4  = 0.739 

q5  = 2.  0 

= -2.  6x  10 

a5  = 0. 440 

r j = 3. 2 

0DO  = 0.  963x10* 4 

P =10.0 

♦ = -105.0 

r2  = “4* 4 

= -0.  101x10" 

0 = -182.5 

r3  a -1.  Ox  10” 5 

u>2  * 0.198x10* 3 

02  = -90.0 
bQ  = 1.368 

r4  = 1. 5xl0"5 

u)3  a *0.  853x10"* 

uk  = 0.106xl0x8 
4 

dQ  = 0. 2784 

bj  = 0.589 

al  = -1.139 

dj  = 0. 1263 

= -15.25 

b2  = 0. 449  xlO”1 

<*2  - 113*0 

d2  = 0.6422x10" 

« = -5.563 

-1 

b3  = 0.468x10 

»3  = 41 . 08 

d3  = 0.3222x10" 

*3  = -1.458 

f 0=  15.57 

eQ  = 0.1149 

Cj  = -19.86 

fj  = 0.6066 

Dj  = 176.5 

ex  = 0.4306x10* 

«2  = 97.99 

f2  = 0.  2784 

*>2  = 17. 32 

e2  = 0. 1186x10" 

«3  = 51.45 

f3  = 0. 2574 

e3  = 68.02 

e3  = 0. 1739x10" 

gQ  * -0.1236 

r\l  = l.  379 

hQ  * 1.  511 

Yx  = 2.221 

gl  - 1. 112 

r\2  = 7.  242 

hx  = 1.325 

V2  = 102.4 

g2  * 0.2338 

n3  = 59.09 

h2  = 0.  3508 

Y3  = 2.  900 

g3  = 0. 2562 

h3«  0.2319 

TABLE  II- 1 
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Both  the  ITS  model  and  the  AFCRL  polar  model  are  based  largely  on 
‘he  vertical  incidence  ionosonde  data  of  IGY  1958.  The  AFCRL  model 
differs  from  the  ITS  model  in  that  the  daily-hourly  observations  of  foF2 
were  used  rather  than  the  monthly  median  data  for  the  set  of  available 
ionosonde  locations.  The  scale  of  significant  polar  ionospheric  phenomena 
both  spatially  and  temporally  is  smaller  than  that  of  the  available  data  and 
the  variability  of  foF2  as  a function  of  the  relevant  parameters  is  so  great 
that  the  data  may  be  insufficient  to  specify  a model  with  the  correct  empir- 
ical behavior.  The  use  of  the  daily  data  allowed  for  a description  of  the 
planetary  magnetic  activity  index,  Kp  dependence  of  foF2  in  the  AFCRL 
model.  This  dependence  was  coupled  to  the  behavior  of  the  important 
ionospheric  boundaries  in  the  polar  region.  The  improved  ITS  model  is 
still  based  upon  median  data  and,  therefore,  does  not  contain  any  Kp 
dependence.  While  the  first  ITS  model  specifies  foF2  on  a global  scale, 
the  AFCRL  model  does  not  contain  a provision  for  predicting  foF2  below 
4fiP  north  geographic  latitude.  This  tends  to  be  a severe  limitation  of  the 
model  since  many  applications  require  knowledge  of  the  ionospheric  param- 
eters over  a radio  wave  path  spanning  different  global  regions. 

Empirical  Basis  for  foF2  Model 

In  order  to  determine  the  behavior  of  foF2,  several  sets  of  ionosonde, 
data  were  used.  Approximately  300, 000  hourly  values  of  foF2  for  the  fifty 
stations  north  of  4£^  latitude  for  IGY  1958  constituted  the  primary  data 
base.  In  addition,  we  used  the  hourly  icnsonde  data  for  each  day  of  March, 
June,  September  and  December  of  1964.  For  some  of  the  study,  we  used 
the  hourly  data  for  the  stations  listed  in  table  II- 2 for  the  years  I960,  1962, 
1966  and  1968. 

The  polar  ionosphere  may  be  described  in  terms  of  several  boundaries 
which  display  certain  regular  characteristics.  The  most  notable  ones  being 
the  trough  and  the  auroral  oval.  .In  an  attempt  to  isolate  these  features  from 


TABLE  II- 2 


Stations  for  Which  foF2  Data  is  Available  for  I960,  1962,  1966,  1968 


Stations 

X (Degrees) 

Xjn  (Degrees) 

0 (Degrees 

Gorky 

156 

56.2 

51.6 

44.3 

Julensruk 

055 

54.6 

62.5 

13.4 

Kiruna 

167 

67.8 

64.3 

20.5 

Leningrad 

160 

60.0 

55.8 

30.3 

Lindau 

050 

51.4 

48.3' 

10.1 

Lycksele 

164 

64.  6 

61.2 

18.8 

Miedzeszyn 

152 

52.  2 

48.2 

21.2 

Moscow 

155 

55.5 

51.0 

37.3 

Murmansk 

168 

68.9 

64. 6 

33.0 

Nurmijarvi 

159 

60.5 

56.  6 

24*  6 

Providenya 

664 

64.4 

60.3 

18.5 

Sakekhard 

266 

66.6 

61.8 

66.  6 

Slough 

051 

51.5 

49.9 

359.4 

Sodankyla 

166 

67.4 

63.4 

26.  6 

Uppsala 

158 

59.8 

56.4 

17.6 

Yakutsk 

462 

62.  0 

56. 4 

129.8 
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the  other  dynamic  processes  occurring  in  the  ionosphere,  we  have  found  it 
useful  to  examine  the  data  for  various  geographic  regions  separately.  We 
have  defined  our  boundaries  by  the  three  parameters:  8^,  the  southern 
boundary  (in-corrected  geomagnetic  latitude)  of  the  auroral  oval,  or  equiva- 
lently the  northern  wall  of  the  trough,  the  northern  boundary  of  the 
auroral  oval  and  0g,  the  southern  boundary  of  the  trough.  These  parameters 
are  given  by  the  following: 

Magnetic  Day  6.  (Xt^^  <18.0  (where  is  the  corrected  geomagnetic 

time) 


0A  = 78.  0°  - 2. 5 Kp  (Degrees) 
0t,  = 82.0  - 1.5  K 

B p 

9S  = 0A-  6.  CP 


Magnetic  Night  ^GM^6’0  0r  tCGM>18‘° 


9A  = 71.0° -E.5K  cos  {jL(tCGM-l)}) 


0B  = 75.0°  + <1.0  - 2.5  cos  {&  (tCGM-l)} 


K 


SS 


Using  these  boundaries,  we  refer  to  the  data  as  belonging  to  one  of  the 
following  regions: 


a)  Xm  < 8^  South  of  the  trough  (ST) 

b)  0g  £ Xm  L 0A  in  the  trough  (T) 
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c)  \ inthe  auroral  oval  (A) 

d)  0j5  < Xm  north  of  the  auroral  oval  (NA) 

where  \m  is  the  corrected  geomagnetic  latitude  of  the  station.  We  also 
treat  the  data  as  occurring  during  the  solar  day  when  y (solar  zenith  angle) 

€.  94. 6°  and  solar  night  y > 94.  6°.  An  indication  of  the  distribution  of  foF2 
observations  can  be  seen  from  the  1958  data.  For  the  daytime  auroral 
region  the  fraction  of  total  observations  for  each  month  varies  between  a 
low  of  2%  during  the  winter  months  and  a peak  of  11%  during  the  summer 
months.  For  the  nighttime  during  the  summer  months,  1%  and  during  the 
winter  months  7%  of  the  data  is  from  the  auroral  oval  region.  For  the 
trough,  this  fraction  of  monthly  data  varies  between  a low  of  2%  during  the 
summer  months  to  14%  during  the  winter  months.  These  small  figures  for 
1958  still  represent  a reasonable  amount  of  data  for  the  study  of  the  behavior 
of  foF2  in  the  two  regions  of  interest.  The  problem  of  having  enough  useful 
data  is  enhanced  during  the  other  years  under  study.  The  distribution  of 
data  for  I960  is  shown  in  the  following  table. 
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TABLE  II -3 


Spatial  Distribution  of  Hourly  foF2  Observations  for  I960 
Number  of  daily-hourly  observations. 


(ST) 

Solar  Dav 
(T) 

' (A); 

<NA) 

(ST) 

Solar  Night 
(T)  (A) 

(NA) 

January 

2900 

250 

0 

0 

4300 

700 

50 

0 

February 

3600 

300 

0 

0 

3500 

600 

50 

0 

March 

5000 

400 

30 

0 

2900 

650 

50 

0 

April 

5000 

800 

100 

0 

1800 

550 

150 

0 

May 

6700 

1200 

150 

0 

1300 

250 

40 

0 

June 

6800 

1400 

100 

0 

1000 

75 

0 

0 

July 

6800 

1200 

125 

0 

1200 

100 

10 

0 

August 

6000 

900 

75 

0 

1700 

325 

50 

0 

September 

5400 

480 

40 

.0 

2700 

500 

60 

0 

October 

4000 

400 

100 

0 

3100 

550 

100 

0 

November 

3100 

200 

30 

0 

3600 

600 

60 

0 

December 

2600 

200 

15 

0 

4100 

750 

90 

0 

It  should  be  noted  that  the  placement  of  each  observation  depends  upon  the 
definition  of  the  boundaries  for  that  hour  and  hence  upon  the  Kp  value.  We 
use  for  our  Kp  data  base  the  three  hour  planetary  index.  We  did  not  attempt 
to  establish  our  boundary  definitions  from  the  data.  Instead,  we  have  con- 
centrated our  efforts  in  trying  to  extract  the  behavior  of  the  amplitude  of 
foF2  in  the  various  regions. 

Corrections  to  AFCRL  Polar  Model 

The  AFCRL  Polar  Model  for  foF2  was  an  attempt  to  parameterize  the 
behavior  of  foF2  in  a limited  geographic  region  in  terms  of  the  spatial  and 
temporal  variables  which  appear  to  be  related  to  the  dynamic  processes 
occurring  in  the  polar  region.  As  can  be  seen  from  the  model  description. 
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they  can  be  divided  into  the  following  categories: 

(1)  Geographic  latitude  and  local  time 

(2)  Geomagnetic  latitude  and  geomagnetic  time 

While  the  diurnal  variation  of  the  model  did  not  predict  the  sunrise  gradients 
with  great  accuracy,  the  largest  deviations  between  the  model  and  the  1958 
data  were  due  to  the  terms  depending  upon  the  magnetic  coordinates.  It  was 
decided  to  redetermine  that  part  of  the  model  for  foF2  which  represents  the 
magnetic  effects.  We  used  as  a starting  model,  foF2m,  that  part  of  the  polar 
model  which  remains  when  the  term 

G (Xm,  m(hp),  D)  = 0 

and  when  the  two  corrections  for  magnetic  day  and  night  are  omitted.  An 
examination  of  the  residuals 


AfoF2  (t)  = foF2™  (t)  - OB S (t) 

at  the  various  ionosonde  locations  showed  the  presence  of  a strong  auroral 
oval  during  the  winter  months.  The  trough  was  not  as  clearly  defined.  The 
residuals  were  then  assumed  to  be  of  the  form: 

AfoF2  = (cq  + Cj  B)  + (c£  + Cj  B + m (c^  + Cg  B))  Xm 


2 

+ (c6  + c?  B + m (cg  + c9  B))  Xm 
+ (CjQ  + cn  B) 

+ (t-to)  (Cj2  ci3  B)  + (ci4  + ci5  ^ ^C16  * c17  ^ ^m 
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+ (c18  + c19 


} 


+ 


(c20  + C21  B)  + 


(c22  + c23  B + m (c24  + c25  3))  *m 


+ (c^^  4 ^27  B ^ ^ (®9A  4 c^q  B))  X 


28  T 29  ""  m 


+ (c30  + c 


31>  Xm) 


where  B = cos  2^B  + ®)  1 ajid  m = m (kp) 


A multiple  linear  regression  procedure  was  used  with  all  the  daily  hourly 
observations  of  foF2  for  the  polar  stations  for  1958.  The  results  of  the 
analysis  were  not  found  to  be  satisfactory.  An  examination  of  the  relative 
regression  coefficients  c^/c^,  Cg/c3,  Cg/cg,  c^/c^,  etc.  indicated  that  no 
dependence  upon  kp  could  be  established  from  this  procedure.  The  kp 
indices  represent  three  hour  global  averages  and  while  they  may  appear  to 
be  correlated  with  the  variation  of  foF2  at  a particular  station  show  little 
correlation  with  foF2  at  sets  of  ionosonde  locatiqps.  The  deficiency  of  the 
regression  technique  is  attributed  to  the  nonuniformity  and  sparseness  of 
the  data  in  the  geographical  region  of  interest.  The  residual,  AfoF2  is 
expected  to  describe  the  variation  of  foF2  in  the  region  of  the  trough  and 
the  auroral  oval  but  as  we  have  already  shown  less  than  10%  ox'  the  data 
are  from  these  regions.  The  regression  procedure  could  only  be  of  value 
if  variable  weighting  factors  could  be  assigned  to  each  of  the  observations. 

Corrections  to  the  ITS  Median  Model 

The  methods  employed  in  the  development  of  the  ITS  median  model  for 
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foF2  tended  to  give  more  weightto  low  and  midlatitude  observations  (ban 
to  those  at  high  latitude.  It  has  been  shown  that  the  re  stilting  numerical 
maps  do  not  represent  the  true  physical  variations  of  the  original  data  in 

(3) 

the  polar  region  . Attempts  have  been  madecto  improve  the  representa- 
tion of  ionospheric  characteristics  in  the  polar  region  using  different 
mapping  coordinates,  but  with  basically  the  same  mapping  functions. 

An  improvement  in  the  representation  of  foF2  in  the  polar  region  was 
obtained  but  the  resulting  model  is  applicable  only  above  30°  north  latitude. 

Since  the  original  global  model  represents  the  physical  behavior  of 
the  ionosphere  quite  well  away  from  the  polar  region,  we  have  chosen  to 
use  this  model  as  a first  approximation  to  the  ionosphere  at  high  latitudes. 
We  have  then  sought  to  find  corrections  to  this  model  to  describe  the  be- 
havior of  the  high  latitude  ionosphere  in  terms  of  the  relevant  spatial  and 
geophysical  variables.  In  an  attempt  to  isolate  the  variation  of  foF2  with 
magnetic  activity  on  a global  scale,  we  first  studied  the  data  available 
during  the  solar  day  (solar  zenith  angle,  x*94. 6°).  For  each  hourly  ob- 
servation of  foF2,  we  have  calculated  the  percent  residual: 


A(t,  0,  X,  s , kp)  «= 


ITS  obs 

foF2  (t,  0,  X,s)  - foF2 

ITS 

foF2  (t,  0,  X,  s) 


were  s is  the  mean  sunspot  number  for  a given  month.  In  figure  II- 1,  we 
show  A as  a function  of  m(kp)  (as  described  in  the  CRL  polar  model)  for 
each  of  the  months  for  the  1958  data.  The  four  areas,  ST,  T,  A,  and  NA, 
of  the  polar  region  are  shown  as  separate  curves.  The  eight  months  from 
March  through  October  show  a similar  linear  dependence  upon  m(kp)  which 
is  independent  of  Xm  for  each  month.  The  winter  months  also  exhibit  the 
same  behavior  except  for  the  region  of  the  auroral  oval  in  which  the  re- 
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RESIDUALS  FOR  1958:  SOLAR  DAY 


Figure  II- lb 


45 


quired  correction  is  much  larger  than  for  the  other  regions.  It  is  also 
important  to  note  that  there  is  very  little  data  available  during  the  solar 
day  for  the  winter  months  so  that  some  of  (these  effects  may  not  be  statis- 
tically significant.  The  data  in  each  ofthesemonths  has  been  separately 
fit  to  the  function: 


A*  = *4  + b|m(kp)  where  i denotes  the  month. 

In  figure  H-2,  we  plot  a^  and  as  a function  of  the  month.  Both  a and  b 
appear  to  have  a six-month  period.  Superimposed  upon  the  experimental 
curves  are  the  two  expressions: 

a(D)  = -0.  155  + 0. 095  cos  (-^^) 

b(D)  = 0.  0425  - 0.  0225  cos(^|2) 

where  D has  been  taken  to  be  the  middle  day  of  each  month.  For  the  low 
solar  activity  year,  19&4  only  four  months  of  data  were  available.  Figure 
H-3  shows  the  variation  of  the  residual,  A>  with  m(kp)  for  the  solar  day  for 
these  four  months.  Here,  we  also  note  only  a very  small  difference  be- 
tween the  separate  areas  within  the  polar  region.  If  we  assume  that  the 
period  and  phase  of  the  a and  b parameters  are  the  same  for  the  both 
years,  we  can  determine  the  amplitudes  for  1964.  The  correction  param- 
eters are: 


and 


a(D)  = -0.  133  + 0.078  cos(~~) 
b(D)  = 0. 0225  - 0. 01 10  cos  (4^2 ) 
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There  is  fairly  good  agreement  between  the  two  years  for  each  of  the  four 
parameters.  The  data  from  each  of  the  other  years  displays  a similar  be- 
havior with  the  value  of  a being  much  more  irregular  than  the  behavior  of 
b.  We  consider  this  to  represent  a global  correction  to  the  ITS  model  for 
both  day  and  night  which  should  be  applied  independent  of  the  polar  region. 
This  given  as: 

foF2(x*94.6°)  = foF2ITS  { 1 - a(D)  - b(D)  m(kp)}-  . 


We  have  examined  the  residuals  for  the  solar/ night  for  each  of  the 
years  as  a function  of  month,  corrected  magnetic  latitude,  corrected 
magnetic  time  and  magnetic  activity.  While  the  residuals  are  clearly  a 
function  of  the  corrected  magnetic  latitude  it  is  difficult  to  extract  any 
regular  variation  with  the  other  parameters.  In  general,  we  have  estab- 
lished the  residual  in  die  area  south  of  the  trough  as  a standard  for  any 
given  month  and  magnetic  activity  and  sought  to  obtain  corrections  in  the 
other  regions  that  are  consistent  with  the  dominant  morophological  features 
of  the  ionosphere.  Our  modeling  of  both  the  trough  and  auroral  oval 
regions  are  strongly  correlated  with  our  assumptions  about  the  location 
and  extent  of  each  of  these  regions  and  their  dependence  upon  magnetic 
activity.  It  has  been  observed  previously  ^ ^ that  the  data  from  ground 
based  ionosondes  during  1958  do  not  show  strong  evidence  of  the  trough. 

We  find  this  to  be  true  even  when  the  dependence  of  the  nighttime  residuals 
on  magnetic  activity  is  al'-o  examined. 

When  the  features  do  occur,  we  expect  to  observe  the  residuals  in  the 
trough  region  to  be  larger  than  those  of  the  region  used  as  a standard. 
Correspondingly,  we  expect  to  observe  residuals  which  are  lower  (ITS 
model  to  small)  in  the  auroral  oval  region  ths^  those  in  standard  region. 

We  illustrate  periods  in  which  these  effects  are  present  with  strong  intensity 
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in  the  following  graphs  of  the  residual  as  a function  of  magnetic  activity 
(m(kp)j,> 

(1)  March  and  June  1964  figure  II- 4.  (X  > 94.  6°) 

(2)  March  and  August  1962  figure  II- 5. 


In  each  of  these  months,  both  features  are  observed.  We  have  found  evi- 
dence for  a nighttime  auroral  oval  correction  in  most  of  the  months  with  it 
being  more  intense  during  the  winter.  The  evidence  for  the  trough  is  much 
more  sporadic  although  it  appears  to  be  quite  weak  during  the  winter.  The 
data  does  not  indicate  a strong  correlation  of  the  intensity  of  the  two  effects 
with  the  three-hour  planetary  index,  k . 

r 

From  the  general  features  of  the  residuals,  we  have  assumed  correc- 
tions to  the  ITS  model  during  the  solar  night  to  be  of  the  form: 

foF2  = foF2(x*94.6°)  £ 1 + &N]  where 
Aq  yjexpxj,  exp  {-  x \/2  }for  XJn^0A 


0A  = 70.  2°  - 1.65  kp  - 5. 1 cos  jj~  _1)]  (degrees) 

XA  ~ 3*  7^  + I*  3 kp  (degrees) 

1 " XA 

^ = T [l  + cos  for  Xm<«A 
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T = Tj  \Jexp  Xj  exp  Xj/2}  exp  (^CGM  ' 3)2/l 2} 


Tj  = 0 for  x s 90° 

°r  6‘0<tCGM<18-0 


T1  = T0  X>94.6°  and  6. 0 < tCGM  < 18.  0 

Ti  = T°  {"■  4 ^ } 90°ix  s 94. 6°  and  6.0<tCGM<18.0 

and  where  tCGM  is  the  corrected  geomagnetic  time  and  \m  is  the  corrected 

geomagnetic  latitude.  The  maximum  amplitudes  have  initially  been  taken 

to  be.  Ty  - -0.  20  and  AQ  = 0.  20.  An  indication  of  the  magnitude  of  the 

correction  in  the  trough  region  is  given  in  table  II- 4 where  the  maximum 

nighttime  values  are  given  for  X = 0.  and  t -in 

m *A  lCGM  " 


TABLE  II-4 

Ratio  of  foF2  / foF2  (ITS) 
NfiP_ 


Day  of  Year 


2 


A 


For  other  values  of  Xm  and  the  ratio  1.  between  the  night  and  day 

values.  The  diurnal  variation  of  foF2  is  shown  in  figure  II- 6.  The  curves 
represent  the  prediction  of  the  ITS  model  without  corrections  and  with  correc- 
ts ons  for  m(kp)  values  of  2,  4,  and  6 for  two  stations  during.  January  1958.  . 

In  addition  to  the  spatial  independent  kp  corrections  applied  to  both  locations, 
the  trough  correction  is  clearly  apparent  in  the  prediction  for  Ottawa. 


Figures  II-7  and  II-8  are  contour  plots  of  the  predictions  of  foF2  for  the 

two  polar  models  for  the  region  north  of  40°  N latitude  for  January  15,  1958 

UT  = 000  hours.  The  division  between  solar  day  and  solar  night  is  explicitly 

shown  (x  =94.6°).  The  second  figure  contains  the  corrections  to  the  ITS 

model  that' are  described  above  for  a k value  of  1.0.  The  ITS  model  is  a 

P 

function  of  universal  time  and  geographic  longitude  whereas  our  polar  region 
corrections  are  also  functions  of  local  time.  The  entire  region  shows  a gen- 
eral increase  in  foF2  due  to  the  k correction.  It  should  be  noted  that  file 

P 

model  is  continuous  in  foF2  and  its  first  derivations  everywhere  on  the  globe. 
We  are  now  in  the  process  of  comparing  the  model  with  additional  foF2  data 
with  an  aim  toward  a better  determination  of  some  of  the  correction  terms. 
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DIURNAL  VARIATION  OF  £oF2:  ITS  MODEL  WITH  POLAR  CORRECTIONS 
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IH.  IMPROVEMENTS  IN  RAY  TRACING  TECHNIQUES 


A Three  Dimensional  Ray  Tracing  Computer  Program  (ARCON  II) 


In  the  past  year,  we  have  continued  our  work  in  the  development  of 
improved  ray  tracing  techniques.  This  work  encompassed  several  major 
modifications  to  the  ARCON  Version  I three  dimensional  ray  tracing 
program.  ^ These  changes  allowed  for  greater  versatility, 

faster  execution  times,  more  efficient  computer  core  allocation,  reduction 
in  programming  required  by  the  user  and  the  consolidation  of  alternative 
subprograms  implementing  various  user  options. 

The  following  pages  describe  the  subprograms  in  which  substantive 
changes  were  made  to  ARCON  Version  I.  No  substantive  changes  were 
made  to  the  subroutines:  TRACE,  HASEL,  RKAM,  BACKUP,  REACH  and 
POLCAR.  The  job  setup  is  described  in  TABLE  HI-1.  Details  of  the  input 
file  are  described  in  TABLE  IH-2  and  the  details  of  the  ionospheric  specifi- 
cation file  are  described  ia. TABLE  HI-3. 
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A.  SUBROUTINE  RAYTRC 


RAYTRC  is  the  main  subroutine  for  the  ARCON  ray  tracing  program. 

It  calls  the  initialization  routine  READW,  and  organizes  calls  to  the  ray 
tracing  subroutine  RANGER  in  one  of  three  sequencing  modes.  These 
modes  are  "scatter, " "homing,"  and  "skip  distance"  sequencing  which 
are  respectively  selected  by  a "1",  "2"  or  "3"  punch  in  column  30  of  the 
option  card  described  in  the  INPUT  file  description. 

"Scatter"  sequencing  is  the  original  sequencing  mode.  Loops  are 
constructed  to  step  by  equal  intervals  through  ranges  in  elevation,  azimuth, 
and  frequency.  Certain  modifications  were  implemented  to  allow  for  desired 
options  in  particular  applications.  First,  die  scan  in  elevation  angles  is 
terminated  as  soon  as  a ray  in  the  set  has  penetrated  the  ionosphere.  This 
feature  may  be  turned  off  by  specifying  a negative  elevation  step  in  the  W 
array,  or  by  punching  character  "1"  in  column  80  of  the  option  card. 

Second,  should  the  azimuth  step  be  zero  subsequent  rays  in  the  elevation 
scan  are  corrected  for  azimuthal  deviation  to  attempt  closer  down  range 
landing  at  a given  azimuth.  This  feature  may  be  turned  off  by  punching 
character  "1"  in  column  75  pf  the  option  card,  or  by  specifying  a non  zero 
azimuth  step.  Third,  an  algorithm  to  compute  the  backscatter  amplitude 
at  each  frequency  was  included.  This  calculation  is  not  accomplished  unless 
at  least  two  elevations  and  two  azimuths  are  traced. 

"Homing"  sequencing  implements  a linear  Newton  iteration  algorithm 
in  elevation  and  azimuth  to  determine  the  elevation  and  azimuth  at  each 
frequency  to  illuminate  a given  down  range  target.  Since  the  algorithm  is 
linear,  it  tends  to  be  unstable  near  the  skip  distance.  For  the  first  fre- 
quency, the  initial  elevation  is  taken  from  the  W array  (element  255)  and 
the  initial  azimuth  is  taken  as  the  azimuth  to  the  target.  The  user  should 
take  care  to  specify  the  initial  elevation  less  than  the  skip  distance  elevation 
and,  if  F mode,  transmission  is  desired,  greater  than  the  E to  F mode 
switching  elevation.  The  initial  elevation  for  the  second  frequency  is  taken 
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as  the  solution  elevation  for  the  first  frequency.  For  subsequent  frequencies, 
the  initial  elevations  are  extrapolated  from  the  solutions  for  the  previous  two 
frequencies.  The  initial  azimuths  for  the  second  and  subsequent  frequencies 
are  taken  as  the  solutions  for  the  last  frequencies.  The  user  specifies  the 
geographic  latitude  and  east  longitude  as  elements  (11)  and  (12),  respectively, 
of  the  W-array.  The  program  computes  the  required  range  and  azimuth. 
Analysis  of  the  homing  algorithm  implemented,  and  discussion  of  its  applica- 
tion in  various  ionospheres  is  contained  in  Section  IV  of  this  report. 

"Skip  Distance"  sequencing  quickly  locates  the  skip  distances  for  the 
sequence  of  frequencies  at  a single  azimuth.  Elevations  are  scanned  back- 
ward from  the  final  elevation  (W(256))  toward  the  initial  elevation  (W(255))  in 
1°  steps  until  the  ionosphere  ceases  to  be  penetrated;  thereafter,  the  elevation 
scan  continues  towards  lower  elevations  in  1/2°  steps  until  the  ground  range 
begins  to  increase.  The  elevation  is  incremented  by  the  elevation  step  (W(257)) 
the  frequency  is  incremented  by  the  frequency  step  (W(260))  and  the  process 
is  repeated  until  either  the  frequency  loop  is  completed  or  an  elevation  angle 
less  than  the  initial  elevation  (W(255))  is  encountered. 
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B.  SUBROUTINE  READW 


Initialization  features  originally  scattered  among  several  different 
subroutines  were  collected  together  in  the  INPUT  file  processing  routine 
READW.  This  program  initializes  labelled  common  blocks  with  data 
statements,  clears  the  INPUT  file  of  the  next  run  identification  card, 
optic,  card,  and  W array  cards.  Non  zero  elements  of  the  W array  are 
listed,  and  elements  of  the  W array  are  distributed  to  other  labelled 
common  blocks.  The  orthogonal  transformation  matrices  among 
geographic,  dipole,  and  transmitter  coordinate  systems  are  computed 
via  calls  to  XFRM.  The  coordinates  of  each  pole  and  receiver  in  each 
system  are  listed.  Elements  of  the  option  list  are  checked  to  construct 
a message,  listed  later,  describing  the  options  for  that  run.  A call  to 
FILEIN  initializes  the  ionosphere.  The  time  required  for  initialization 
tasks  is  computed  and  written  to  the  output  file. 
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C.  COMPLEX  FUNCTION  XFRM  (X,  R) 


This  is  a general  application  subprogram  to  compute  Hie  orthogonal' 
transformation  matrix  between  various  coordinate  systems.  Alternate 
entry  points  GTOD  and  DTOG  compute  the  transformed  spherical  coordinates. 

The  argument  R (nine  element  array)  is  the  orthogonal  transformation 
matrix,  while  X (declared  complex)  contains  colatitude  and  east  longitude 
in  the  real  and  imaginary  parts  respectively.  Similarly,  the  function 
return  contains  colatitude  and  east  longitude. 

Entry  point  XFRM  computes  Hie  transformation  matrix  R to  a coordinate 
system  with  pole  through  point  X in  Hie  old  system.  The  prime  meridian  of 
the  new  system  is  the  great  circle  containing  both  poles.  The  function 
return  is  the  new  pole  X. 

Entry  point  GTOD  returns  the  coordinates  in  the  new  system  corres- 
ponding to  coordinates  X in  the  old,  while  DTOG  returns  old  system 
coordinates  corresponding  to  new  system  coordinates  X. 
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D.  SUBROUTINE  FILEIN  (N,  M,  PFA) 


FILEIN  processes  the  ionosphere  file  described  in  Table  IH-1.  The 
arguments  specify  which  subfile  to  read  from,  how  many  records  to  list,  > 
and  where  to  store  the  plasma  frequency  array: 

N < O Rewind  file,  position  to  (-N)^1  subfile. 

N > O No  rewind,  position  to  subfile. 

N = O No  file  action,  continue  with  current  ionosphere. 

M < O Single  profile  option.  Read  and  list  (-M)**1  profile, 
skip  to  end  of  subfile,  set  up  ionosphere  domains  to 
reflect  universal  spherically  symmetric  ionosphere. 

M O Read  entire  file,  list  fii  st  M+l  profiles. 

PFA  First  word  address  of  plasma  frequency  array. 

The  routine  first  positions  the  ionosphere  file  to  the  desired  subfile, 
and  reads  the  ionosphere  parameter  record.  The  upper  and  lower  limits 
of  the  domain  of  the  ionosphere  in  kilometers  of  altitude,  degrees  of  latitude, 
and  degrees  of  longitude  are  computed  and  stored  in  the  first  six  elements 
of  common  block/ CRDLMTS/.  The  ionosphere  parameters  are  printed. 

The  core  required  to  store  the  array  is  computed  and  a call  to  NRFL 
establishes  the  field  length  at  the  required  element  of  PFA.  The  ionosphere 
is  then  read  in  profile  by  profile,  and  listed  as  requested. 

Note  that  FILEIN  does  not  require  positive  definite  grid  spacing s,  a 
feature  which  removes  constraints  on  the  user  previously  imposed. 

Note  also  that  the  user  is  relieved  of  the  burden  to  allocate  core 
memory;  function  NRFL  (described  below)  dynamically  increases  or 
decreases  field  length  to  reflect  requirements. 
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E.  FUNCTION  NRFL(A) 


This  is  a FORTRAN  callable  CDC6600  COMPASS  subprogram  that 
allows  the  user  to  increase  or  decrease  core  memory  allocated  to  the 
job  at  any  time  during  execution,  subject,  of  course,  to  the  maximum 
core  limit  specified  on  the  job  card. 

The  argument  A must  be  a simple  variable  name  or  array  element 
in  blank  common.  The  subprogram  establishes  the  field  length  of  the  job 
at  the  address  defined  by  the  argument,  writes  a message  to  the  dayfile, 
FL-nnnnnn,  where  nnnnnn  is  the  address  defined  by  A,  and  returns  via 
the  function  call  the  address  of  A. 
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F.  SUBROUTINE  PRINTR  (NWHY) 

Changes  resulting  in  substantial  core  memory  savings  were  imple- 
mented in  the  basic  output  routine.  These  consisted  primarily  in  removing 
the  initialization  features  to  the  initialization  subprogram.  The  output  line 
format,  however,  was  not  changed.  The  argument  NWHY  still  remains  as 
a display  code  menomic  indicating  the  ray  states  triggering  the  call  to 
PRINTR.  Each  output  line  contains  the  following  information: 

Item  Format  Description 


1 

IXE8. 1 

relative  error  in  group  speed  squared 

2 

IXA8 

NWHY 

3 

F8.2 

altitude  in  kilometers 

4 

F9.2 

ground  range  in  kilometers 

5 

F7.  2 

angle  in  degrees  between  ray  and  field 

6 

F7.2 

magnetic  colatitude 

7 

F7.2 

magnetic  east  longitude 

8 

F7.2 

azimuthal  deviation 

9 

F7.2 

local  azimuthal  deviation 

10 

F7.2 

elevation  of  field  point  from  transmitter 

11 

F7.2 

local  elevation  of  ray  path 

12 

F7.2 

real  part  of  polarization 

13 

F7.2 

imaginary  part  of  polarization 

14 

F8.2 

group  path  length 

15 

F8.2 

phase  path  length 

The  first  call  to  PRINTR  for  each  ray  traced  will  print  column 
headings.  PRINTR  will  recompute  the  magnetic  field  for  purposes  of 
computing  item  5 if  the  output  magnetic  field  (field  7 of  the  option  card) 
specifies  a more  accurate  magnetic  field  than  included  in  the  refractive 
index  calculation. 
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G.  FUNCTION  HANGER  (FREQ,  EL,  AZ) 

This  function  subprogram  was  inserted  as  a buffer  between  the  ray 
selection  routine  RAYTRACE  and  the  ray  tracing  routine  TRACE.  The 
three  arguments  are,  respectively,  ray  frequency  in  megahertz,  takeoff 
elevation  in  radians,  and  takeoff  azimuth  in  radians.  The  function  return 
is  the  ground  range  in  kilometers  (taken  from  the  PRINTR  output  list)  upon 
return  from  subroutine  TRACE. 

The  program  performs  some  initialization  for.  each  ray  traced.  It 
writes  one  line  on  the  output  file  specifying  the  frequency  in  megahertz, 
azimuth  in  degrees,  and  elevation  in  degrees.  It  computes  the  initial 
conditions  for  the  ray,  and  aborts  if  the  ray  is  in  an  evanescent  region. 

It  calls  in  the  fundamental  ray  tracing  routine  TRACE,  and  prints  the  time 
required  to  execute  each  ray. 


H.  SUBROUTINE  ELECTX 


This  subroutine  computes  the  parameter  X and  its  derivatives,  where 
X is  the  squared  ratio  of  ionosphere  plasma  frequency  to  working  frequency. 

The  plasma  frequency  and  its  derivatives  e computed  by  calling  a 
general  second  order  interpolation  algorithm  ii  .aree  variables,  FOXYZ. 
Initialization  of  the  common  blocks  required,  and  input  of  the  plasma  fre- 
quency array  is  accomplished  in  SUBROUTINE  FILEIN. 

The  program  must  convert  the  field  point  variables  from  geocentric 
distance  in  kilometers  to  altitude  in  kilometers,  colatitude  in  radians  to 
latitude  in  degrees,  and  east  longitude  in  radians  to  east  longitude  in  degrees. 
Furthermore,  an  initial  domain  check  of  longitude  is  necessary,  adding  or 
subtracting  360°  if  out  of  range,  before  calling  FOXYZ.  This  removes 
from  the  user  the  burden  of  unnecessary  rigor  of  selecting  longitude  branch 
domain. 

Should  the  user  desire  to  use  an  analytic  electron  density  model,  he 
need  only  replace  this  routine  and  SUBROUTINE  FILEIN  (Nl,  N2,  PFA). 

The  user  may  use  the  latter  routine  for  initialization  purposes;  it  is  called 
once  for  each  W-array  specified  on  the  INPUT  file. 
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I.  FUNCTION  FOXYZ(X,  Y,  Z,  F,  IX,  IY,  IZ) 


i 


i 


! 


FOXYZ  implements  a general,  second  order  interpolation  algorithm 
within  a tabulated  function  of  three  variables.  The  elements  of  the 
argument  list  are: 

X,  Y,  Z - coordinates  of  the  field  point  at  which  the  interpolated 
value  of  the  function  is  desired 

F - first  word  address  of  tabulated  array  of  function 

values 

IX,  IY,  IZ  - array  dimensions 
The  function  return  is  the  interpolated  function  value. 


Prior  to  calling  this  routine,  elements  of  certain  labelled  common 
blocks,  /CRDLMTS/and/IONOSPH/,  must  be  defined  in  addition  to  filling 
in  the  array  F.  These  elements  are  listed  below: 


/ CR DLMTS / XMIN ( 3) , XMAX(3),  MINMAX(3),  IFLG,  LN(2) 
/IONOSPH/NAME,  LN(3),  XO(3),  DX(3),  DM(5) 

XMIN  lower  limits  of  function  definition 
XMAX  upper  limits  of  function  definition 
MINMAX  defined  by  FOXYZ,  see  below 
IFLG  defined  by  FOXYZ,  see  below 
not  used 


LN 

NAME 

LN 

XO 

DX 

DM 


not  used 

array  dimensions,  must  agree  with  IX, IY,  IZ  of  argument  lis 
coordinates  of  first  element  of  F 
uniform  grid  spacings 
not  used 


In  the  above  arrays,  indices  1,  2,  3 correspond  respectively  to 
coordinates  X,  Y,  Z;  thus,  the  proper  domain  of  varh-,ble  X is  given  by 


XMIN(l)  < X ^ XMAX(l). 
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The  program  first  checks  the  field  point  variables  against  their 
respective  proper  domains.  If  any  variable  is  less  than  its  minimum 
value,  the  corresponding  element  of  MINMAX  is  set  to  -1,  if  greater 
than  the  maximum,  MINMAX  is  set  to  +1;  otherwise,  MINMAX  is  set 
to  zero.  If  any  variable  is  outside  its  proper  domain,  IF LG  is  set  to 
+ 1,  otherwise,  IFLG  is  set  to  zero. 

If  variable  X is  outside  its  proper  domain,  a value  of  0 is  returned 
via  the  function  call.  However,  if  Y and/or  Z are  outside  their  domains, 
second  order  interpolation  proceeds  only  in  X,  at  the  nearest  profile. 

The  program  stores  local  variables  in  labelled  common  block /SCRATCH/, 
64  words  in  length,  certain  elements  of  which  may  be  of  interest  and  use 
to  the  user.  These  elements  are  listed  below. 

/SCRATCH/SC(64) 

element  definition  on  return  from  FOXYZ 


12,13,14  dF/dx,  dF/Sy,  3F/3z 
15,16,17  a2F/3x2,  32F/3y2,  32F/aZ2 

18,19,20  a2F/ayaz,  a2F/azax,  a2F/axay 


J.  SUBROUTINE  RINDEX 


Thii?  program  was  modified  to  include  all  cases.  The  program  calls 
subprograms  ELECTX,  DIPY,  and  COLFRZ  to  compute  the  electron 
density,  magnetic  fields,  and  electron  collision  frequencies  respectively. 
Each  contribution  is  value  checked,  and  if  zero  complicated  computations 
are  bypassed  to  execute  simpler  algorithms. 

Core  storage  conservation  was  accomplished  by  (1)  sharing  local 
variable  storage  with  other  routines  in  labelled  common  block  /SCRATCH/, 

(2)  using  DO  loops  for  contributions  by  the  magnetic  field  components, 

(3)  eliminating  non  essential  complex  variable  calculations,  and  (4)  elimi- 
nating the  complex  attribute  from  variables  of  which  only  the  real  part  is 
ultimately  used.  Many  of  these  steps  also  resulted  in  reduction  of  execution 
time. 


K.  FUNCTION  DIPY  (YEAR,  R,  THETA,  PHI,  NORDER) 

This  subprogram  computes  a spherical  harmonic  approximation  to  the 
earth*  s magnetic : field.  The  arguments  are: 

YEAR  - yea?  to  which  to  update  the  harmonic  coefficients 

R - geocentric  distance  in  kilometers 

THETA  - magnetic  dipole  colatitude  in  radians 

PHI  - magnetic  dipole  east  longitude  in  radians 

NORDER  - order  of  approximation  from  O (no  field)  to  8 (maximum) 

The  function  return  is  the  magnetic  potential. 


The  components  and  their  derivatives  of  the  magnetic  field  are  returned 
via  labelled  common  block/DIPYOUT/,  eighteen  words  in  length: 

word  definition 

1 magnetic  potential 

2 field  magnitude  in  gammas 

3, 4, 5 derivatives  of  field  magnitude 

6-9  radial  component  and  derivatives 

10-13  southward  component  and  derivatives 

14-17  eastward  component  and  derivatives 

1 8 dip  angle 


TABLE  in-1:  JOB  SETUP  FOR  RAYTRACE  PROGRAM 

JOB  CARD 

ATTACH,  TAPE4.  Cards  to  attach  file  of  ionospheres 
ATTACH,  MODE  LX,  MODELX3141,  ID=GIBBS,  MR=1. 

LIBRARY,  MODE  LX. 

FTN. 

LGO. 

7/8/9 

PROGRAM  MAIN  (INPUT,  OUTPUT,  TAPE4) 

CALL  RAYTRC 
END 

7/8/9 

IDENTIFICATION  CARD 

OPTION  CARD 

W ARRAY  CARDS 

• 

• 

BLANK  CARD  - causes  execution 

(Repeat  as  often  as  desired.  Terminate  run  with  W-array  element 
negative) 

6/7/8/9 
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TABLE  IQ-2:  INPUT  FILE  EXPECTED  BY  READ  W 


CARD  FORMAT 


1 

2 

3 toN 
N+J 


FORMAT  DESCRIPTION 

8*}c°  descriptive  text  to  be  listed  on  OUTPUT 

1615  option  card 

13,  E14.7,4Il  W array  cards 

Repeat  from  card  1 as  often  as  desired 


OPTION  CARD  DESCRIPTION 
FIELD  FUNCTION 


USED  IN 


^0  produces  output  to  TAPE  1 PRINTR 

= 0 inhibits  output  to  OUTPUT  PRINTR 

selects  ionosphere  from  TAPE  4 FILEIN 
selects  profiles  to  list  FILEIN 

selects  magnetic  field  option  DIPY  & RINDEX 

selects  sequencing  mode  RAYTRC 

selects  output  magnetic  field  PRINTR 

reserved  for  collision  frequency  option 
not  used 

^ 0 inhibits  azimuth  correction  RAYTRC 
= 0 inhibits  elevation  scan  RAYTRC 


WARRAY  CARD  DESCRIPTION  („e  reference  2) 
FIELD  VALUE  FUNCTION 


1-400 

0 

<0,  >400 

REAL 

7*0 

7*0 

7*0 

7*0 


W array  element 

end  of  W- array,  start  ray  trace 

terminate  job 

W array  value 

Convert  value  from  degrees  to  radians 
Convert  value  from  kilometers  to  radians 
Convert  value  from  nautical  miles  to  kilometers 
Convert  value  from  feet  to  kilometers 
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TABLE  III- 3:  TAPE4-ianosphere  file 

STRUCTURE  OF  IONOSPHERE  FILE  expected  by  FILEIN 

FILE  STRUCTURE 


RECORD 

DESCRIPTION 

LENGTH 

- 

1 

IONOSPHERE  GRID  PARAMETERS  (IGP) 

15 

PLASMA  FREQUENCY  PROFILES  (PFP) 

U6 

2 

first  latitude,  first  longitude 

t 

M+l 

last  latitude,  first  longitude 

M+2 

first  latitude,  second  longitude 

M*N+1 

last  latitude,  last  longitude 

M*N+2 

EOF  (may  be  absent) 

0 

M*N+3 

IGP-  second  file 

15 

etc. 


i 


RECORD  STRUCTURE 


WORD 

DESCRIPTION 

TYPE 

DIMENSION 

IGP  1 

Name  of  Ionosphere 

Display  Code 

_ 

2 

Number  of  altitudes 

(L) 

integer 

- 

3 

Number  of  latitudes 

(M) 

integer 

- 

4 

Number  of  longitudes  (N) 

integer 

- 

5 

First  altitude 

(HO) 

real 

kilometers 

6 

First  latitude 

real 

degrees 

7 

First  longitude 

real 

degrees 

8 

altitude  step 

(DH) 

real 

kilometers 

9 

latitude  step 

real 

degrees 

10 

longitude  step 

real 

degrees 

11 

Date  Time  Group 

Display  Code 

- 

12 

F10.7 

read. 

13 

SUNSPOT  number 

real 

14 

AP 

real 

15 

KP 

real 

PFP  1 

Profile  length 

(L) 

integer 

- 

2 

First  altitude 

(HO) 

real 

kilometers 

3 

Profile  latitude 

real 

degrees 

4 

Profile  longitude 

real 

degrees 

5 

last  altitude 

(HM) 

real 

kilometers 

6 

altitude  step 

(DH) 

real 

kilometers 

7 

Plasma  Frequency  at  HO 

real 

MHZ 

U6 

Plasma  Frequency  at  HM 

real 

MHZ 
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IV.  ANALYSIS  OF  A HOMING  TECHNIQUE  IN  THREE  DIMENSIONAL 
IONOSPHERIC  RAY  TRACING 


A.  Development  of  Homing  Algorithm 

We  desire  to  determine  the  takeoff  elevation  and  azimuth  of  a ray 
illuminating  a given  target  from  a given  transmitter  at  a given  frequency. 
Let  the  position  of  the  target  be  specified  by  the  range  R«j>  and  azimuth  Arj. 
from  the  transmitter.  The  relevant  results  from  a ray  trace  at  elevation 
E and  azimuth  A are  the  range  R(E,  A)  and  downrange  azimuth  Aj^  (E,  A). 
The  homing  solution,  E^j  and  A^,  is  found  by  inverting  the  equations. 


R (ER,  -A-jj)  = Rrp, 


(1) 


ar  (eh'  ah)  = at* 

In  the  neighborhood  of  Eq,  Aq,  R (E,  A)  and  AR  (E,  A)  may  be  approx 
imated  by 


R (E,  A)  = Rq  + Re  (E-Eq)  + Ra  (A-Aq) 

Ar  (E,  A)  = Arq  + AE  (E-E0)  + Aa  (A-Aq) 


(2) 


where  Rq  = R (Eg,  Ag),  A^g  = Aj^(Eg,Ag),  Rg  and  R^  are  the  derivatives 
of  R at  Eg,  Ag,  and  Ag,  AA  are  the  derivatives  of  Ar. 

Thus  the  inversion  of  equations  reduces  to  salving  for  Epj  and  Ajj  the 
linear  simultaneous  equations, 
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(4) 


% (eh  ' =0)  + RA  <AH  ‘ A0>  = rt  " R0  = AR 
AE  <%  ” E0>  + AA  <AH  " A0^  = AT  ” AR0  = 

Thus, 

EH  + E0  + AE# 
ah  = A0  + ^ 

where, 

AE  = (AR.  aa  - AAr-  ra)  / (re-  aa  - ae-  Ra) 

AA  = (AAr**Re  - AR*  Ae)  / (re‘  aa  - Ae*  Ra)  . 

Starting  from  E , A as  a good  first  approximation,  iteration  of 
equations  (4)  converges  rapidly  to  the  solution. 

In  principle,  computation  of  the  derivatives  appearing  in  equations  (5) 
requires  tracing  three  rays  at  each  iteration.  However,  considerable 
simplification  is  possible. 

The  ranges  of  a family  of  ray  tracings  are  primarily  dependent  on 
elevation  and  relatively  insensitive  to  azimuth  changes.  Furthermore, 
the  downrange  azimuthal  deviation  6 is  slowly  varying  and  relatively  in- 
sensitive to  changes  in  initial  conditions.  Thus,  we  may  set  AE  = RA  = 0, 
and  A.  = 1;  equations  (2)  simplify  to, 


R (E,  A)  = Rq  + Re  (E  - Eq) 
A (E,  A)  = A - 6 


(7) 
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Equations  (4)  for  successive  approximations  to'ihe  homing  solution  reduce 
to, 


Ejj  " E0  + ~ V ^ RE 

Ajj  = Aij  + 6 


(8) 


The  remaining  derivative  of  R is  given  by  the  finite  difference  formula. 


RE  = (R(E1>  - r(e2))/  (Ej  - e2) 


(9) 


The  final  algorithm  implemented  in  the  homing  mode  of  ray  tracing 
generates  a sequence  of  elevations  and  azimuths  E^,  given  by 

E(0)  = E0  (10) 

E(1)  = E*0)  + AE0 

E<n+1)  = E<n)  + (E(n)  - E(n_10  ’ (rt  ’ r(”0  / CH(n)  - R(n"10 
A<°>  =AT 

A(n+1)  = AT  + 6 (n) 

where  Eg  and  AEq  are  specified  by  the  user,  and  5^  are  the  range  and 
azimuthal  deviation  of  the  nth  iteration.  The  user  specifies  geographic 
latitude  and  east  longitude  of  the  target  receiver  from  which  the  program 
computes  A^  and  R^.. 

Iteration  of  equations  (10)  terminates  when  the  range  error  |Rij«  - Rl  '| 
is  less  than  one  kilometer  or  when  n = 10. 


-79- 


B.  Ray  Tracing  in  Rush-Miller  Ionospheres 

The  Rush- Miller  Electron  Density  Model  with  ITS  model  £oF2  and  HznF2 

was  used  to  generate  ionospheres  of  plasma  frequency  profiles  over  the 

domain  36°  N to  4°  N magnetic  and  14°  W to  55°  E magnetic.  The  solar  cycle 

was  taken  at  mean  conditions  (solar  flux  100  and  sun  spot  number  50),  and 

ll  ll 

equinox  (March  21)  ionospheres  were  constructed  for  0 and  6 Eastern 
Standard  Time.  These  two  ionospheres  will  be  referred  to  as  the  Midnight 
and  Sunrise  ionospheres  hereafter. 

We  traced  rays  in  both  ionospheres  from  a transmitter  located  at  24.  5°N 
and  81.5°W  geographic  at  an  azimuth  of  108°  True  and  a frequency  of  6 MHZ. 
This  frequency  is  about  20%  larger  than  the  mean  F2  critical  frequency  along 
the  path.  We  used  a null  magnetic  field,  and  integrated  the  ray  tracing  equa- 
tions with  8 kilometer  steps. 

Figure  IV- 1 displays  plasma  frequency  profiles  for  several  ranges  along 
azimuth  108°  True  for  the  sunrise  ionosphere,  and  shows  the  development  of 
the  ionospheric  E layer  from  shorthy  before  to  about  an  half  hour  after  sun- 
rise. 

Figure  IV- 2 is  a presentation  of  several  ray  traces  across  the  sunrise 
line.  The  altitude  scale  is  four  times  larger  than  the  range  scale.  These 
are  crude  ray  traces  in  that  only  the  critical,  points  are  plotted  (enter  and 
exit  ionosphere,  reflection  and  ground  return)  connected  by  straight  lines. 

Figure  IV- 3 is  a range  versus  elevation  presentation  of  ray  traces  in 
the  midnight  and  sunrise  ionospheres. 

Even  though  the  sunrise  line  is  positioned  from  the  transmitter  to  present 
the  most  difficulty,  figures  IV-2  and  IV-3  reflect  considerable  regularity.  Note 
from  figure  IV-,2  that  the  switch  over  from  E to  F mode  transmission  (3°  to  4°) 
occurs  while  the  rays  are  entering  the  ionosphere  on  the  sunlit  side. 

The  major  anomaly  in  ray  tracing  across  the  sunrise  line  occurs  for 
rays  (14°  to  17°)  entering  the  ionosphere  near  the  high  altitude  knee  of  the 


-80- 


FI&IJRE  EC-l 

PLASMA  FREQUENCY  PROFILES 
AT  SEVERAL  RANGES,  FROM  TRANSMITTER 
ALONG  AZIMUTH  lOS’TRUE 


ma 


I OOOYjtA 


F XCi’Arh  3Z-^ 


F LEV/.  7 . r // 


i 


JS.-H 

RA-V^r  vs,  Fu-vw.  ••  • 

PASS'.  •:  . f'J  4 J,  ,•'//,•;  JHU'.’FtA'.G  */*,•£ 

EteQHe*/cY  t>  h*± 

X »■  faRAIWCr  4PP«C/*4K  £/^TH 


* 


/; 

< t 


/ x 


* 

/i 

/ i 

-i  i 

i 
! 


%V>£.v»moiy 

53  n>°  )5° 


— j 


84 


ionosphere  bottom.  The  high  horizontal  gradients  at  this  point  tend  to  turn 
the  ray  upwards,  and  effectively  create  a hole  in  the  ionosphere  through 
which  the  rays  penetrate  into  space.  This  anomaly,  of  course,  reflects  the 
structure  of  the  Rush-Miller  ionospheric  model  which  may  or  may  not  be  a 
fully  adequate  representation  of  the  actual  conditions  that  prevail  near  sun- 
rise and  sunset. 

Although  the  ray  trace  program  appears  to  be  reliable  for  tracing  rays 
across  the  sunrise  line,  severe  difficulties  associated  with  extreme  trans- 
verse gradients  are  encountered  when  tracing  rays  along  the  sunrise. 

Figure  IV-4  represents  range  elevation  curves  from  a transmitter  located 

0 

1 east  of  sunrise  transmitting  close  to  the  sunrise  line. 

At  azimuth  180°  True  transmission  is  entirely  in  the  daylight  side.  The 
range  elevation  curve  is  relatively  smooth.  Rays  are  deflected  westward  by 
3°  to  4°  in  F-mode  transmission. 

At  azimuth  200°  True,  transmission  is  entirely  in  the  night  ionosphere 
up  to  11°  elevation.  Beyond  11°  the  ray  enters  the  day  ionosphere,  exits  the 
ionosphere  across  the  sunrise  surface  and  reenters  the  night  ionosphere. 

The  range  elevation  curve  (Figure  IV-4)  is  close  to  that  at  180°  demonstrating 
the  dominance  of  the  F2  layer  in  determining  transmission  in  F mode. 

At  azimuth  190°  True,  transmission  is  along  the  sunrise  surface. 

Figure  IV-4  shows  that  the  ray  tracing  program  does  not  yield  stable  results. 
No  improvement  in  smoothing  this  curve  resulted  on  reducing  the  integration 
step  size  from  8 kilometers  to  1 kilometer.  Clearly,  this  instability  would 
preclude  meaningful  results  in  an  attempt  to  compute  a homing  solution. 


C.  Homing  in  Rush-Miller  Ionospheres 

The  homing  mode  of  the  ray  trace  program  was  used  to  compute  homing 
solutions  from  a transmitter  at  24.5°N,  8L5°W  geographic  to  a target 

receiver  at  20.  0°N,  69.  0°W  geographic  at  u frequency  of  6 MHZ.  This 
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corresponds  to  a ground  range  of  1379  km  at  an  azimuth  of  108.  8°  True. 

The  initial  elevation  was  selected  from  the  range  elevation  curves  of  figure 
IV- 3. 

The  results  of  successive  rays  in  approaching  the  solution  are  tabulated 
in  tables  IV- 1 and  IV -2  for  homing  in  the  midnight  and  sunrise  ionospheres 
respectively.  The  appropriate  points  are  also  plotted  on  the  range  elevation 
curves  of  figure  IV- 3. 

In  file  midnight  ionosphere  convergence  to  the  homing  solution  was 
rapid.  The  transverse  gradients  in  .he  ionosphere  are  small  resulting  in 
lees  than  4 kilometers  in  azimuthal  deviation  at  the  neighborhood  of  the 
target. 

Homing  across  the  sunrise  surface  tended  to  converge  more  slowly  the 
down  range  azimuthal  deviations  corresponded  to  corrections  at  the  target 
of  about  20  kilometers  reflecting  the  higher  transverse  gradients  in  the 
ionosphere. 

As  expected,  an  attempt  to  home  nearly  parallel  to  the  sunrise  surface 
was  unsuccessful  ■**  the  inst«uultty  of  the  range  elevation  curve. 


TABLE  IV-!:  HOMING  IN  MIDNIGHT  IONOSPHERE 


Frequency:  6 MHZ  No  magnetic  field 

Receiver  Range:  1379  Km  Azimuth:  108.  8°True 


Ray  Nr 

Takeoff 

Azimuth 

(deg) 

Takeoff 

Elevation 

(cleg) 

Azimuthal 

Deviation 

(deg) 

Range 
Error 
(K  m ) 

1 

108. 80 

11. 00 

+ 0.  12 

-426. 93 

2 

108.  91 

12.  96 

+ 0.  13 

-266. 27 

3 

108. 93 

16.  20 

+ 0.  12 

-51. 87 

4 

108. 92 

16.  98 

+ 0.  20 

+ 3.60 

5 

109.  00 

16.  93 

fO.  20 

-0.  75 

TABLE  IV- 2:  HOMING  ACROSS  SUNRISE  SURFACE 


Frequency:  6 MHZ 
Receiver  Range:  1379  Km 

No  magnetic  field 
Azimuth:  108.8°  True 

Ray  Nr 

Takeoff 

Azimuth 

(deg) 

Takeoff 

Elevation 

(deg) 

Azimuthal 

Deviation 

(deg) 

Range 

Error 

(Km) 

J 

108. 80 

11.  00 

-1.23 

+98. 14 

2 

107.  ‘37 

10.  55 

-0. 

337.01 

3 

107. 84 

10.  28 

-0.  96 

3 21 .44 

4 

107. 84 

9.  90 

-0.  86 

-8.  74 

5 

107. 93 

10.  01 

-0.  91 

+ 1.62 

6 

107. 89 

10.  00 

-0.  88 

- 3.  36 

7 

107. 92 

10.  01 

0.  89 

3 0.  70 

